#include "cem/constants.h"
#include "cem/electromagnetics_system.h"
#include "cem/parameter_reader.h"
#include "cem/post_processor.h"

#include <libmesh/dof_map.h>
#include <libmesh/elem.h>
#include <libmesh/enum_order.h>
#include <libmesh/equation_systems.h>
#include <libmesh/fe.h>
#include <libmesh/fe_interface.h>
#include <libmesh/linear_implicit_system.h>
#include <libmesh/mesh_base.h>
#include <libmesh/numeric_vector.h>
#include <libmesh/point_locator_tree.h>
#include <libmesh/quadrature_gauss.h>
#include <libmesh/vector_value.h>
#include <libmesh/vtk_io.h>

#include <cmath>
#include <complex>
#include <fstream>
#include <iostream>
#include <limits>
#include <map>
#include <vector>

namespace cem
{

// 函数: write_basic_log
// 作用: 打印系统名与解向量的基本统计，若存在 NaN/Inf 则列出前若干个异常条目
void
PostProcessor::write_basic_log(const ElectromagneticsSystem & system) const
{
    // 并行信息与系统标识
    auto & es = const_cast<libMesh::EquationSystems &>(system.equation_systems());
    const int rank = es.get_mesh().comm().rank();

    std::cout << "[POST][rank " << rank << "] enter write_basic_log for system='" << system.name()
              << "'" << std::endl;

    // 获取线性系统与解向量原始指针（unique_ptr -> raw pointer）
    const auto & lin = system.linear_system();
    const auto * sol_ptr = lin.solution.get();
    if (!sol_ptr)
    {
        std::cout << "[POST][rank " << rank << "] solution vector is null" << std::endl;
        return;
    }

    // 打印尺寸信息（可能涉及全局/本地）
    std::cout << "[POST][rank " << rank
              << "] solution sizes: global=" << static_cast<unsigned long long>(sol_ptr->size())
              << ", local=" << static_cast<unsigned long long>(sol_ptr->local_size()) << std::endl;

    // 基本范数（对复数向量，libMesh 定义为幅值相关的范数）
    const double l2_norm_value = sol_ptr->l2_norm();
    const double linfty_norm_value = sol_ptr->linfty_norm();

    // 本地化到主机内存以便进一步统计（复数幅值的min/max、NaN/Inf等）
    std::vector<libMesh::Number> local_values;
    try
    {
        sol_ptr->localize(local_values);
    }
    catch (const std::exception &)
    {
        local_values.clear();
    }
    catch (...)
    {
        local_values.clear();
    }

    double min_abs_value = std::numeric_limits<double>::infinity();
    double max_abs_value = 0.0;
    std::size_t nan_count = 0, inf_count = 0;
    const std::size_t print_limit = 10;
    std::size_t printed_count = 0;

    if (!local_values.empty())
    {
        for (std::size_t i = 0; i < local_values.size(); ++i)
        {
            const auto & z = local_values[i];
            const double abs_v = std::abs(z);

            if (std::isfinite(abs_v))
            {
                if (abs_v < min_abs_value)
                    min_abs_value = abs_v;
                if (abs_v > max_abs_value)
                    max_abs_value = abs_v;
            }
            else
            {
                if (!std::isfinite(std::real(z)) || !std::isfinite(std::imag(z)))
                {
                    if (std::isnan(std::real(z)) || std::isnan(std::imag(z)))
                        ++nan_count;
                    if (!std::isfinite(std::real(z)) || !std::isfinite(std::imag(z)))
                        ++inf_count;

                    if (printed_count < print_limit)
                    {
                        std::cout << "[POST] non-finite entry at index=" << i << ", value=" << z
                                  << std::endl;
                        ++printed_count;
                    }
                }
            }
        }
    }

    bool has_bad_entry = !(std::isfinite(l2_norm_value) && std::isfinite(linfty_norm_value)) ||
                         nan_count || inf_count;

    std::cout << "[POST] solution stats => size="
              << static_cast<unsigned long long>(sol_ptr->size()) << ", l2_norm=" << l2_norm_value
              << ", linfty_norm=" << linfty_norm_value << ", min_abs=" << min_abs_value
              << ", max_abs=" << max_abs_value << std::endl;

    if (nan_count || inf_count)
    {
        std::cout << "[POST] WARNING: solution non-finite entries => NaN=" << nan_count
                  << ", Inf=" << inf_count << std::endl;
    }

    if (!has_bad_entry)
    {
        std::cout << "[POST] solution appears finite." << std::endl;
    }
}

// 函数: write_vtk_snapshot
// 作用: 使用 libMesh::VTKIO 将当前方程系统（包含网格与解）写出为 VTK/并行VTK 文件（由后缀决定）。
// 说明: - 通过 VTKIO::write_equation_systems 写出 .pvtu（并行）或 .vtu/.vtk（串行）
//       - 写出前调用 EquationSystems::update()，确保 current_local_solution 与派生场同步
void
PostProcessor::write_vtk_snapshot(const ElectromagneticsSystem & system,
                                  const std::string & filename) const
{
    try
    {
        // 写出前先做一次解摘要检查
        write_basic_log(system);

        auto & es = const_cast<libMesh::EquationSystems &>(system.equation_systems());
        es.update();

        // 基于底层网格构建 VTK 写出器
        libMesh::VTKIO vtk_io(es.get_mesh());
        vtk_io.write_equation_systems(filename, es);

        std::cout << "[POST] VTK written -> " << filename << std::endl;
    }
    catch (const std::exception & ex)
    {
        std::cerr << "[POST] VTK write failed: " << ex.what() << std::endl;
    }
}

// 函数: write_solution_csv
// 作用: 将解向量以 CSV 文本格式写出（index,value_real,value_imag）。仅 rank 0 写文件。
void
PostProcessor::write_solution_csv(const ElectromagneticsSystem & system,
                                  const std::string & filename) const
{
    // 功能: 直接导出单复变量解向量的实部/虚部，格式: index,u_r,u_i（仅 rank 0 写文件）。
    auto & es = const_cast<libMesh::EquationSystems &>(system.equation_systems());
    const int rank = es.get_mesh().comm().rank();

    const auto & lin = system.linear_system();
    const auto * sol_ptr = lin.solution.get();
    if (!sol_ptr)
    {
        if (rank == 0)
            std::cerr << "[POST] solution vector is null, skip writing '" << filename << "'"
                      << std::endl;
        return;
    }

    using number_t = libMesh::Number;
    std::vector<number_t> values;
    try
    {
        sol_ptr->localize_to_one(values, 0);
    }
    catch (...)
    {
        if (rank == 0)
            std::cerr << "[POST] localize_to_one failed, cannot write CSV" << std::endl;
        return;
    }

    if (rank == 0)
    {
        std::ofstream ofs(filename);
        if (!ofs.is_open())
        {
            std::cerr << "[POST] cannot open file '" << filename << "' for writing" << std::endl;
            return;
        }
        ofs << "index,u_r,u_i\n";
        for (std::size_t i = 0; i < values.size(); ++i)
        {
            ofs << i << "," << std::real(values[i]) << "," << std::imag(values[i]) << "\n";
        }
        ofs.close();
        std::cout << "[POST] solution CSV written -> " << filename << std::endl;
    }
}

// 作用: 按 FE 解在采样点直接计算 MT 阻抗 Zxy=Ex/Hy 与 Zyx=Ey/Hx，
//       并换算视电阻率 rho_app=|Z|^2/(mu0*omega) 与相位（度）。
// 说明: 当前实现仅在 rank 0 上执行，并仅对 rank 0 所拥有的本地活动单元上的点进行插值；
//       若点不在本地单元内，则输出 NaN 以提示用户（可后续扩展为并行聚合）。
void
PostProcessor::write_mt_csv(const ElectromagneticsSystem & system,
                            const ParameterReader & params,
                            const std::string & filename) const
{
    // 并行支持：所有进程共同参与点定位与数值计算，仅在 rank 0 聚合输出到文件
    auto & es = const_cast<libMesh::EquationSystems &>(system.equation_systems());
    const libMesh::MeshBase & mesh = es.get_mesh();
    const auto & comm = mesh.comm();
    const int rank = comm.rank();

    const auto & lin = system.linear_system();
    const libMesh::DofMap & dof_map = lin.get_dof_map();
    const unsigned int var = system.var_index();

    // FE 构建：与系统变量一致的 Nedelec 元；设置一个占位的 Gauss 积分（自定义点会覆盖）
    const unsigned int dim = mesh.mesh_dimension();
    const libMesh::FEType fe_type = lin.variable_type(var);
    std::unique_ptr<libMesh::FEVectorBase> fe = libMesh::FEVectorBase::build(dim, fe_type);
    libMesh::QGauss qrule(dim, libMesh::SECOND);
    fe->attach_quadrature_rule(&qrule);

    // 取形函数与 curl 访问器（注意: 仅在调用 reinit 后有效）
    const auto & phi_v_all = fe->get_phi();
    const auto & curl_v_all = fe->get_curl_phi();

    // 点定位器（物理点 -> 单元）
    std::unique_ptr<libMesh::PointLocatorBase> locator =
        std::make_unique<libMesh::PointLocatorTree>(mesh);
    // locator->init();

    // 常数与材料映射
    const double omega = system.omega();
    const double mu0 = cem::constants::vacuum_permeability;
    const auto & mu_map = system.material_manager().area_mu_map(); // 值为相对磁导率 mu_r

    // 采样点坐标
    cem::CollectPointValues pts = params.compute_collect_point_values();
    if (pts.x_coords.empty())
        pts.x_coords = {0.0};
    if (pts.y_coords.empty())
        pts.y_coords = {0.0};
    if (pts.z_coords.empty())
        pts.z_coords = {0.0};

    // 访问解向量（用于按局部自由度抽取）
    using number_t = libMesh::Number;
    const auto * sol_ptr = lin.solution.get();
    if (!sol_ptr)
    {
        if (rank == 0)
            std::cerr << "[POST] solution vector is null, skip MT export" << std::endl;
        return;
    }

    // 仅在 rank 0 打开输出文件
    std::ofstream ofs;
    if (rank == 0)
    {
        ofs.open(filename);
        if (!ofs.is_open())
        {
            std::cerr << "[POST] cannot open file '" << filename << "' for writing" << std::endl;
            return;
        }
        ofs << "num,x0_loc,y0_loc,z0_loc,xy_ris,yx_ris,alph_xy,alph_yx\n";
    }

    // 遍历采样点
    std::size_t idx = 0;
    const double tiny = 1e-20;

    std::vector<libMesh::dof_id_type> dofs;
    std::vector<number_t> vals;

    for (double x : pts.x_coords)
        for (double y : pts.y_coords)
            for (double z : pts.z_coords)
            {
                const libMesh::Point p_phys(x, y, z);
                const libMesh::Elem * elem = (*locator)(p_phys);

                double rho_xy = std::numeric_limits<double>::quiet_NaN();
                double rho_yx = std::numeric_limits<double>::quiet_NaN();
                double phase_xy_deg = std::numeric_limits<double>::quiet_NaN();
                double phase_yx_deg = std::numeric_limits<double>::quiet_NaN();

                bool ok = false;
                if (elem && !elem->is_remote())
                {
                    // 求参考单元坐标
                    libMesh::Point pref =
                        libMesh::FEInterface::inverse_map(dim, fe_type, elem, p_phys);

                    // 用自定义点重入 FE
                    std::vector<libMesh::Point> ref_pts(1, pref);
                    fe->reinit(elem, &ref_pts);

                    // 抽取该单元的单变量自由度与解值(复数)
                    dof_map.dof_indices(elem, dofs, var);
                    const unsigned int n = dofs.size();
                    if (n > 0 && phi_v_all.size() == n && curl_v_all.size() == n)
                    {
                        vals.assign(n, number_t());
                        sol_ptr->get(dofs, vals);

                        // 取 qp=0 的形函数值与 curl 形函数
                        const unsigned int qp = 0;
                        libMesh::RealVectorValue e_r(0., 0., 0.), e_i(0., 0., 0.);
                        libMesh::RealVectorValue curl_e_r(0., 0., 0.), curl_e_i(0., 0., 0.);

                        for (unsigned int i = 0; i < n; ++i)
                        {
                            const auto & phi_i = phi_v_all[i][qp];
                            const auto & curl_i = curl_v_all[i][qp];

                            const double ur = std::real(vals[i]);
                            const double ui = std::imag(vals[i]);

                            e_r.add_scaled(phi_i, ur);
                            e_i.add_scaled(phi_i, ui);
                            curl_e_r.add_scaled(curl_i, ur);
                            curl_e_i.add_scaled(curl_i, ui);
                        }

                        // 由 Faraday 定律: curl(E) = -i * omega * mu * H
                        // 则 H 的实/虚部可由 curl(E) 的虚/实部线性组合得到（仅实标量缩放）
                        const int sub_id = static_cast<int>(elem->subdomain_id());
                        const double mu_r = (mu_map.count(sub_id) ? mu_map.at(sub_id) : 1.0);
                        const double mu = mu0 * mu_r;

                        double inv_om_mu = 0.0;
                        if (omega > 0.0 && std::isfinite(mu) && std::fabs(mu) > tiny)
                            inv_om_mu = 1.0 / (omega * mu);
                        libMesh::RealVectorValue h_r = curl_e_i;
                        h_r *= (-inv_om_mu);
                        libMesh::RealVectorValue h_i = curl_e_r;
                        h_i *= (inv_om_mu);

                        const std::complex<double> ex(e_r(0), e_i(0));
                        const std::complex<double> ey(e_r(1), e_i(1));
                        const std::complex<double> hx(h_r(0), h_i(0));
                        const std::complex<double> hy(h_r(1), h_i(1));

                        const auto make_rho_phase = [&](const std::complex<double> & Z)
                        {
                            // 计算视电阻率
                            const double rho = std::norm(Z) / (mu0 * omega);
                            double phase =
                                -std::atan2(Z.imag(), Z.real()) * 180.0 / cem::constants::pi;
                            return std::pair<double, double>(rho, phase);
                        };

                        if (std::abs(hy) > tiny)
                        {
                            const std::complex<double> zxy = ex / hy;
                            auto rp = make_rho_phase(zxy);
                            rho_xy = rp.first;
                            phase_xy_deg = rp.second;
                        }
                        if (std::abs(hx) > tiny)
                        {
                            const std::complex<double> zyx = -ey / hx;
                            auto rp = make_rho_phase(zyx);
                            rho_yx = rp.first;
                            phase_yx_deg = rp.second;
                        }
                        ok = true;
                    }
                }

                // 跨进程归约：拥有该点的本地进程提供数值，rank 0 负责写出
                int valid_local = ok ? 1 : 0;
                int valid_sum = valid_local;
                comm.sum(valid_sum);

                double rho_xy_sum = ok ? rho_xy : 0.0;
                double rho_yx_sum = ok ? rho_yx : 0.0;
                double phase_xy_sum = ok ? phase_xy_deg : 0.0;
                double phase_yx_sum = ok ? phase_yx_deg : 0.0;
                comm.sum(rho_xy_sum);
                comm.sum(rho_yx_sum);
                comm.sum(phase_xy_sum);
                comm.sum(phase_yx_sum);

                if (rank == 0)
                {
                    if (valid_sum == 0)
                    {
                        std::cerr << "[POST] point (" << x << "," << y << "," << z
                                  << ") not covered by any local element." << std::endl;
                    }

                    const double out_rho_xy = (valid_sum > 0)
                                                  ? (rho_xy_sum / valid_sum)
                                                  : std::numeric_limits<double>::quiet_NaN();
                    const double out_rho_yx = (valid_sum > 0)
                                                  ? (rho_yx_sum / valid_sum)
                                                  : std::numeric_limits<double>::quiet_NaN();
                    const double out_ph_xy = (valid_sum > 0)
                                                 ? (phase_xy_sum / valid_sum)
                                                 : std::numeric_limits<double>::quiet_NaN();
                    const double out_ph_yx = (valid_sum > 0)
                                                 ? (phase_yx_sum / valid_sum)
                                                 : std::numeric_limits<double>::quiet_NaN();

                    ofs << idx++ << "," << x << "," << y << "," << z << "," << out_rho_xy << ","
                        << out_rho_yx << "," << out_ph_xy << "," << out_ph_yx << "\n";
                }
            }

    if (rank == 0)
    {
        ofs.close();
        std::cout << "[POST] MT CSV written -> " << filename << std::endl;
    }
}

void
PostProcessor::write_mt_csv_dual_from_host(const ElectromagneticsSystem & system,
                                           const std::vector<std::complex<double>> & sol_a,
                                           const std::vector<std::complex<double>> & sol_b,
                                           const ParameterReader & params,
                                           const std::string & filename) const
{
    // 并行支持：所有进程共同参与定位与归约，仅在 rank 0 负责写文件
    auto & es = const_cast<libMesh::EquationSystems &>(system.equation_systems());
    const libMesh::MeshBase & mesh = es.get_mesh();
    const auto & comm = mesh.comm();
    const int rank = comm.rank();

    const auto & lin = system.linear_system();
    const unsigned int var = system.var_index();
    const libMesh::DofMap & dof_map = lin.get_dof_map();

    const std::size_t ndofs_global = lin.solution ? lin.solution->size() : 0;

    // 0) 校验输入尺寸，仅在 rank 0 判定，然后广播结果
    int inputs_ok = 1;
    if (rank == 0)
    {
        if (ndofs_global == 0 || sol_a.size() != ndofs_global || sol_b.size() != ndofs_global)
            inputs_ok = 0;
    }
    comm.broadcast(inputs_ok, 0);
    if (!inputs_ok)
    {
        if (rank == 0)
        {
            std::cerr << "[POST] dual-from-host invalid input array sizes (expect " << ndofs_global
                      << ")" << std::endl;
        }
        return;
    }

    // 1) 将 host 上的复数解广播到所有进程（以 real/imag 打包的方式）
    std::vector<double> sol_a_pack(2 * ndofs_global, 0.0), sol_b_pack(2 * ndofs_global, 0.0);
    if (rank == 0)
    {
        for (std::size_t i = 0; i < ndofs_global; ++i)
        {
            sol_a_pack[2 * i + 0] = sol_a[i].real();
            sol_a_pack[2 * i + 1] = sol_a[i].imag();
            sol_b_pack[2 * i + 0] = sol_b[i].real();
            sol_b_pack[2 * i + 1] = sol_b[i].imag();
        }
    }
    comm.broadcast(sol_a_pack, 0);
    comm.broadcast(sol_b_pack, 0);

    std::vector<std::complex<double>> sol_a_all(ndofs_global), sol_b_all(ndofs_global);
    for (std::size_t i = 0; i < ndofs_global; ++i)
    {
        sol_a_all[i] = {sol_a_pack[2 * i + 0], sol_a_pack[2 * i + 1]};
        sol_b_all[i] = {sol_b_pack[2 * i + 0], sol_b_pack[2 * i + 1]};
    }

    const unsigned int dim = mesh.mesh_dimension();
    const libMesh::FEType fe_type = lin.variable_type(var);
    std::unique_ptr<libMesh::FEVectorBase> fe = libMesh::FEVectorBase::build(dim, fe_type);
    libMesh::QGauss qrule(dim, libMesh::SECOND);
    fe->attach_quadrature_rule(&qrule);

    const auto & phi_v_all = fe->get_phi();
    const auto & curl_v_all = fe->get_curl_phi();

    std::unique_ptr<libMesh::PointLocatorBase> locator =
        std::make_unique<libMesh::PointLocatorTree>(mesh);
    // locator->init();

    const double omega = system.omega();
    const double mu0 = cem::constants::vacuum_permeability;
    const auto & mu_map = system.material_manager().area_mu_map();

    cem::CollectPointValues pts = params.compute_collect_point_values();
    if (pts.x_coords.empty())
        pts.x_coords = {0.0};
    if (pts.y_coords.empty())
        pts.y_coords = {0.0};
    if (pts.z_coords.empty())
        pts.z_coords = {0.0};

    std::ofstream ofs;
    if (rank == 0)
    {
        ofs.open(filename);
        if (!ofs.is_open())
        {
            std::cerr << "[POST] cannot open file '" << filename << "' for writing" << std::endl;
            return;
        }
        ofs << "num,x0_loc,y0_loc,z0_loc,xy_ris,yx_ris,alph_xy,alph_yx\n";
    }

    const double tiny = 1e-20;
    std::size_t idx = 0;

    std::vector<libMesh::dof_id_type> dofs;

    for (double x : pts.x_coords)
        for (double y : pts.y_coords)
            for (double z : pts.z_coords)
            {
                const libMesh::Point p_phys(x, y, z);
                const libMesh::Elem * elem = (*locator)(p_phys);

                double rho_xy = std::numeric_limits<double>::quiet_NaN();
                double rho_yx = std::numeric_limits<double>::quiet_NaN();
                double phase_xy_deg = std::numeric_limits<double>::quiet_NaN();
                double phase_yx_deg = std::numeric_limits<double>::quiet_NaN();

                bool ok = false;
                if (elem && !elem->is_remote())
                {
                    libMesh::Point pref =
                        libMesh::FEInterface::inverse_map(dim, fe_type, elem, p_phys);
                    std::vector<libMesh::Point> ref_pts(1, pref);
                    fe->reinit(elem, &ref_pts);

                    dof_map.dof_indices(elem, dofs, var);
                    const unsigned int n = dofs.size();
                    if (n > 0 && phi_v_all.size() == n && curl_v_all.size() == n)
                    {
                        const unsigned int qp = 0;
                        libMesh::RealVectorValue e_a_r(0., 0., 0.), e_a_i(0., 0., 0.);
                        libMesh::RealVectorValue c_a_r(0., 0., 0.), c_a_i(0., 0., 0.);
                        libMesh::RealVectorValue e_b_r(0., 0., 0.), e_b_i(0., 0., 0.);
                        libMesh::RealVectorValue c_b_r(0., 0., 0.), c_b_i(0., 0., 0.);
                        for (unsigned int i = 0; i < n; ++i)
                        {
                            const auto & phi_i = phi_v_all[i][qp];
                            const auto & curl_i = curl_v_all[i][qp];
                            const auto gi = static_cast<std::size_t>(dofs[i]);
                            const std::complex<double> za = sol_a_all[gi];
                            const std::complex<double> zb = sol_b_all[gi];
                            const double ur_a = za.real();
                            const double ui_a = za.imag();
                            const double ur_b = zb.real();
                            const double ui_b = zb.imag();
                            e_a_r.add_scaled(phi_i, ur_a);
                            e_a_i.add_scaled(phi_i, ui_a);
                            c_a_r.add_scaled(curl_i, ur_a);
                            c_a_i.add_scaled(curl_i, ui_a);
                            e_b_r.add_scaled(phi_i, ur_b);
                            e_b_i.add_scaled(phi_i, ui_b);
                            c_b_r.add_scaled(curl_i, ur_b);
                            c_b_i.add_scaled(curl_i, ui_b);
                        }

                        const int sub_id = static_cast<int>(elem->subdomain_id());
                        const double mu_r = (mu_map.count(sub_id) ? mu_map.at(sub_id) : 1.0);
                        const double mu = mu0 * mu_r;

                        // 由 Faraday 定律: curl(E) = -i * omega * mu * H
                        // 则 H 的实/虚部可由 curl(E) 的虚/实部线性组合得到（仅实标量缩放）
                        double inv_om_mu = 0.0;
                        if (omega > 0.0 && std::isfinite(mu) && std::fabs(mu) > tiny)
                            inv_om_mu = 1.0 / (omega * mu);
                        libMesh::RealVectorValue h_a_r = c_a_i;
                        h_a_r *= (-inv_om_mu);
                        libMesh::RealVectorValue h_a_i = c_a_r;
                        h_a_i *= (inv_om_mu);
                        libMesh::RealVectorValue h_b_r = c_b_i;
                        h_b_r *= (-inv_om_mu);
                        libMesh::RealVectorValue h_b_i = c_b_r;
                        h_b_i *= (inv_om_mu);

                        // 从两套极化组装 2x2 系统：E = Z * H，Z = E * H^{-1}
                        const std::complex<double> ex_a(e_a_r(0), e_a_i(0));
                        const std::complex<double> ey_a(e_a_r(1), e_a_i(1));
                        const std::complex<double> hx_a(h_a_r(0), h_a_i(0));
                        const std::complex<double> hy_a(h_a_r(1), h_a_i(1));

                        const std::complex<double> ex_b(e_b_r(0), e_b_i(0));
                        const std::complex<double> ey_b(e_b_r(1), e_b_i(1));
                        const std::complex<double> hx_b(h_b_r(0), h_b_i(0));
                        const std::complex<double> hy_b(h_b_r(1), h_b_i(1));

                        // H 矩阵: [[hx_a, hx_b],[hy_a, hy_b]]，计算其行列式
                        const std::complex<double> detH = hx_a * hy_b - hx_b * hy_a;
                        if (std::abs(detH) > tiny)
                        {
                            // H^{-1} = 1/det * [[ hy_b, -hx_b], [ -hy_a, hx_a ]]
                            const std::complex<double> invH00 = hy_b / detH;
                            const std::complex<double> invH01 = -hx_b / detH;
                            const std::complex<double> invH10 = -hy_a / detH;
                            const std::complex<double> invH11 = hx_a / detH;

                            // Z = E * H^{-1}，其中 E = [[ex_a, ex_b],[ey_a, ey_b]]
                            const std::complex<double> zxx = ex_a * invH00 + ex_b * invH10;
                            const std::complex<double> zxy = ex_a * invH01 + ex_b * invH11;
                            const std::complex<double> zyx = ey_a * invH00 + ey_b * invH10;
                            const std::complex<double> zyy = ey_a * invH01 + ey_b * invH11;

                            const auto make_rho_phase = [&](const std::complex<double> & Z)
                            {
                                const double rho = std::norm(Z) / (mu0 * omega);
                                double phase =
                                    -std::atan2(Z.imag(), Z.real()) * 180.0 / cem::constants::pi;
                                return std::pair<double, double>(rho, phase);
                            };

                            auto rp_xy = make_rho_phase(zxy);
                            // 历史约定：Zyx 取 -Ey/Hx（见单极化实现），等价于在 2×2 解耦后对 Zyx
                            // 取负号
                            auto rp_yx = make_rho_phase(-zyx);
                            rho_xy = rp_xy.first;
                            phase_xy_deg = rp_xy.second;
                            rho_yx = rp_yx.first;
                            phase_yx_deg = rp_yx.second;
                            ok = true;
                        }
                        else
                        {
                            ok = false;
                        }
                    }
                }

                // 2) 对每个采样点进行跨进程归约（拥有该点的本地进程提供数值）
                int valid_local = ok ? 1 : 0;
                int valid_sum = valid_local;
                comm.sum(valid_sum);

                double rho_xy_sum = ok ? rho_xy : 0.0;
                double rho_yx_sum = ok ? rho_yx : 0.0;
                double phase_xy_sum = ok ? phase_xy_deg : 0.0;
                double phase_yx_sum = ok ? phase_yx_deg : 0.0;
                comm.sum(rho_xy_sum);
                comm.sum(rho_yx_sum);
                comm.sum(phase_xy_sum);
                comm.sum(phase_yx_sum);

                if (rank == 0)
                {
                    if (valid_sum == 0)
                    {
                        std::cerr << "[POST] point (" << x << "," << y << "," << z
                                  << ") not covered by any local element." << std::endl;
                    }

                    const double out_rho_xy = (valid_sum > 0)
                                                  ? (rho_xy_sum / valid_sum)
                                                  : std::numeric_limits<double>::quiet_NaN();
                    const double out_rho_yx = (valid_sum > 0)
                                                  ? (rho_yx_sum / valid_sum)
                                                  : std::numeric_limits<double>::quiet_NaN();
                    const double out_ph_xy = (valid_sum > 0)
                                                 ? (phase_xy_sum / valid_sum)
                                                 : std::numeric_limits<double>::quiet_NaN();
                    const double out_ph_yx = (valid_sum > 0)
                                                 ? (phase_yx_sum / valid_sum)
                                                 : std::numeric_limits<double>::quiet_NaN();

                    ofs << idx++ << "," << x << "," << y << "," << z << "," << out_rho_xy << ","
                        << out_rho_yx << "," << out_ph_xy << "," << out_ph_yx << "\n";
                }
            }

    if (rank == 0)
    {
        ofs.close();
        std::cout << "[POST] MT CSV written -> " << filename << std::endl;
    }
}

// 函数: write_real_part_debug
// 作用: 复现 legacy/solver.c#L122-131 的调试输出，写出前 501 个（或更少）非零条目对应的
//       两个实部组成项：matrix_a_r 与 matrix_a_r2*omega*omega。
// 说明: 在新架构中我们不直接持有压缩存储数组，故按 RR 块（实部变量）重建：
//       - matrix_a_r    = ∫_Ω (curl ϕ_i)^T (μ^{-1}) (curl ϕ_j) dΩ
//       - matrix_a_r2   = -∫_Ω (ε ϕ_j)·ϕ_i dΩ   （注意是负号；写出时乘以 ω²）
//       若区域未给出各向异性张量，则按各向同性简化（与装配实现一致）。
void
PostProcessor::write_real_part_debug(const ElectromagneticsSystem & system,
                                     const std::string & filename) const
{
    auto & es = const_cast<libMesh::EquationSystems &>(system.equation_systems());
    const int rank = es.get_mesh().comm().rank();
    if (rank != 0)
        return;

    const auto & lin = system.linear_system();
    const libMesh::DofMap & dof_map = lin.get_dof_map();

    const libMesh::MeshBase & mesh = es.get_mesh();
    const unsigned int dim = mesh.mesh_dimension();

    // FE 准备（与装配一致）
    const unsigned int var = system.var_index();
    const libMesh::FEType fe_type = lin.variable_type(var);
    std::unique_ptr<libMesh::FEVectorBase> fe = libMesh::FEVectorBase::build(dim, fe_type);
    libMesh::QGauss qrule(dim, libMesh::SECOND);
    fe->attach_quadrature_rule(&qrule);

    const std::vector<libMesh::Real> & JxW = fe->get_JxW();
    const std::vector<std::vector<libMesh::RealVectorValue>> & curl = fe->get_curl_phi();
    const std::vector<std::vector<libMesh::RealVectorValue>> & phi_v = fe->get_phi();

    // 材料映射
    const auto & epsl_map = system.material_manager().area_epsl_map();
    const auto & mu_map = system.material_manager().area_mu_map();

    // ω、常数
    const double omega = system.omega();
    const double mu0 = cem::constants::vacuum_permeability;
    const double eps0 = cem::constants::vacuum_permittivity;

    // 使用有序 map 聚合 (row,col) -> (kr, mr2)
    std::map<std::pair<libMesh::dof_id_type, libMesh::dof_id_type>, std::pair<double, double>> acc;

    std::vector<libMesh::dof_id_type> dofs;

    for (auto it = mesh.active_local_elements_begin(); it != mesh.active_local_elements_end(); ++it)
    {
        const libMesh::Elem * elem = *it;

        // 该单元的自由度（仅实部变量）
        dof_map.dof_indices(elem, dofs, var);
        const unsigned int n = dofs.size();
        if (n == 0)
            continue;

        // 绑定 FE 到当前单元
        fe->reinit(elem);

        // 区域标识与材料（各向同性简化）
        const int sub_id = static_cast<int>(elem->subdomain_id());
        const double mu_scalar = (mu_map.count(sub_id) ? mu_map.at(sub_id) : 0.0); // 磁化率 χ
        const double mu_r_eff = mu_scalar;
        const double inv_mu_scalar =
            (std::fabs(mu_r_eff) > 0.0 ? 1.0 / (mu0 * mu_r_eff) : 1.0 / mu0);

        const double epsl_scalar = (epsl_map.count(sub_id) ? epsl_map.at(sub_id) : 0.0); // (ε_r-1)
        const double epsl_eff = (epsl_scalar)*eps0;

        // 积分点循环
        const unsigned int n_qp = qrule.n_points();
        for (unsigned int qp = 0; qp < n_qp; ++qp)
        {
            const double w = JxW[qp];
            for (unsigned int i = 0; i < n; ++i)
            {
                for (unsigned int j = 0; j < n; ++j)
                {
                    const double curl_term = inv_mu_scalar * (curl[i][qp] * curl[j][qp]) * w;
                    const double mass_eps_core = epsl_eff * (phi_v[i][qp] * phi_v[j][qp]) * w;
                    const double k_val = curl_term;
                    const double r2_val = -mass_eps_core; // 注意负号

                    auto key = std::make_pair(dofs[i], dofs[j]);
                    auto it_acc = acc.find(key);
                    if (it_acc == acc.end())
                        acc.emplace(key, std::make_pair(k_val, r2_val));
                    else
                    {
                        it_acc->second.first += k_val;
                        it_acc->second.second += r2_val;
                    }
                }
            }
        }
    }

    // 仅 rank 0 输出前 501 条
    std::ofstream ofs(filename);
    if (!ofs.is_open())
    {
        std::cerr << "[POST] cannot open file '" << filename << "' for writing" << std::endl;
        return;
    }

    std::size_t i = 0;
    for (const auto & kv : acc)
    {
        if (i > 500)
            break;
        const double a_r = kv.second.first;                          // matrix_a_r
        const double a_r2_omega2 = kv.second.second * omega * omega; // matrix_a_r2 * ω²
        ofs << "i= " << i << ", matrix_a_r[" << i << "]= " << a_r << ", matrix_a_r2[" << i
            << "]*omega*omega= " << a_r2_omega2 << "\n";
        ++i;
    }
    ofs.close();
    std::cout << "[POST] Real part debug written -> " << filename << " (entries=" << i << ")"
              << std::endl;
}

} // namespace cem